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Part I 


Introduction 
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From the first airplanes steered by handles, wheels, and pedals to today’s 
advanced aircraft, there has been a century of revolutionary inventions, all of 
them contributing to flight quality. The stability and controllability of aircraft 
as they appear to a pilot are called flying or handling qualities. Many years after 
the first airplanes flew, flying qualities were identified and ranked from desirable 
to unsatisfactory. Later on engineers developed design methods to satisfy these 
practical criteria. 

CONDUIT, which stands for Con trol Designer’s Unified Interface, is a mod- 
ern software package that provides a methodology for optimization of flight con- 
trol systems in order to improve the flying qualities. 

CONDUIT is dependent on an the optimization engine called CONSOL- 
OPTCAD (C-0)[1]. C-0 performs multicriterion parametric optimization. It 
was developed at the Institute for Systems Research (ISR) at the University of 
Maryland at College Park. C-0 was successfully tested on a variety of control 
problems before it was tested on the design of a rotorcraft flight control system 
for the UH-60A helicopter in hover based on the ADOCS controller structure 
[3]. The optimization-based computational system, C-O, requires a particular 
control system description as a MATLAB file and possesses the ability to modify 
the vector of design parameters in an attempt to satisfy performance objectives 
and constraints specified by the designer, in a C-type file. 

After the first optimization attempts on the UH-60A control system, an early 
interface system, named GIFCORCODE (Graphical Interface for CONSOL- 
OPTCAD for Rotorcraft Controller Design) [7] was created. This interface 
allowed the designer to push buttons in order to command the optimization 
engine, instead of typing keywords in the C-0 language. Posing the problem 


2 



still required code writing. 

After that a big step followed, transforming GIFCORCODE into CONDUIT, 
[2], CONDUIT eliminates the code writing, by providing (a) a SIMULINK 
window to design the block diagram, (b) a large library from which to choose 
the handling qualities specifications, (c) a “specmaker” facility to create new 
performance specifications as well as (d) a complete methodology to analyze and 
test the design. Using CONDUIT any control system problem related to aircraft 
can be posed just by using the mouse to click on menus and drag objects from 
libraries. Furthermore, the program provides the designer, after each iteration 
of C-O, with a very efficient picture of all performance criteria plus supporting 
plots for each individual criterion, a plot showing the evolution of the design 
parameters, and a specifications-evolution plot. To benefit from hardware speed 
and advanced visualization characteristics, CONDUIT was developed on the 
Silicon Graphics IRIX operating system. 

CONDUIT has been used to analyze and design: 

• a flight control system for the RASCAL UH-60A helicopter in hover having 
the ADOCS control structure [3] 

• a flight control system for the RASCAL UH-60A in forward flight at 80 
knots [4] 

• a longitudinal control system for the X-29A Forward-Swept Wing Fly-by- 
Wire Demonstrator 

• a roll channel control system for the XV-15 Tiltrotor aircraft in forward 
flight mode 
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• a longitudinal/vertical control system for the Kaman SH-2 rotorcraft in 
hover 

• a lateral/directional control system for the Boeing KC-135 Jet Transport 

• a longitudinal control system for the Grumman F-14 aircraft 

• an antenna control system 

In March 1998 CONDUIT was released to US aircraft and helicopter compa- 
nies, where it is being used in the design of a number of practical aircraft control 
systems. 
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Part II 


Background 



Computer-Aided Control System Design A reasonable way to design a 
control system for a complex large system is by using an optimization technique 
requiring two phases. In the first phase the designer formulates the problem. 
Such a formulation normally consists of an appropriate system structure, a set 
of objectives, and reasonable initial values for the adjustable design parameters. 
During the second phase, an optimization package is used to obtain the best 
values of the design parameters. The second phase is interactive, allowing the 
designer to explore the design tradeoffs. 

Aircraft and helicopter control systems are complex and designing such sys- 
tems usually requires balancing competing objectives. Moreover the process of 
optimizing the performance of the control system for aircraft and helicopters has 
to be repeated tens (or even hundreds) of times corresponding to the changes in 
the system that occur at different flight conditions. 

CONDUIT performs computer-aided control system design using this opti- 
mization technique and it uses C-0[1] as the optimization engine. It is also true 
that many hypothetical designs may be considered — both for the aircraft and 
the controller. 

Specifications for Helicopters and Aircraft Control systems are present in 
almost all modern aircraft and helicopters in order to 1) enhance aircraft mission 
capability, 2) improve handling qualities, and 3) decrease pilot workload. 

Control law design can only be performed satisfactorily if a set of design 
requirements or performance criteria is available. In the case of control systems 
for piloted aircraft generally applicable quantitative design criteria are difficult 
to obtain because the ultimate test of handling qualities is the pilot’s judge- 
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ment. A pilot’s opinion of the flying qualities of an aircraft is influenced by 
the ergonomic design of the cockpit controls, the visibility from the cockpit, 
the weather conditions, the mission requirements, and physical and emotional 
factors. The variability introduced by all these factors can be reduced by aver- 
aging test results over many flights and pilots. A systematic approach to flying 
qualities evaluation is available through pilot opinion rating scales such as the 
Cooper-Harper scale [19]. Once a rating scale has been established it is possible 
to begin correlating the pilot opinion rating with the properties of the aircraft 
dynamic model, and hence derive some analytical specifications that will guar- 
antee good flying qualities. 

The acceptability of flying qualities is quantified in terms of “Levels” that 
are defined for each specific mission task, using the Cooper-Harper scale. The 
Cooper-Harper scale consists of pilot ratings from 1 to 10, where 1 indicates 
“excellent, highly desirable” and 10 means “major deficiencies and lost control 
during some portion of the required operation”. Based on this scale “Levels” are 
defined as follows: Level 1 indicates pilot ratings between 1 and 3 1/2, Level 2 
between 3 1/2 and 6 1/2, and Level 3 between 6 1/2 and 8 1/2. Each Level value 
represents a minimum condition necessary to meet that “Level” of acceptability. 
Level 1 indicates that aircraft characteristics are “good, may be some negligible 
deficiencies” [19], Level 2 shows “moderately objectionable deficiencies” [19], 
and Level 3 shows “major deficiencies” [19] . An aircraft that exhibits Level 
1 flying qualities throughout its flight envelope is fully satisfactory. Some real 
aircraft exhibit some Level 2 flying qualities. Any Level 3 flying qualities are 
unsatisfactory. 
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Chapter 1 


CONDUIT 


The objective of the software package CONDUIT is to act as a designer’s assis- 
tant/ associate for the analysis, testing, and optimization of flight control systems 
in order to improve the flying qualities. In aircraft design, the design specifica- 
tions (and constraints) are often competing and CONDUIT helps the designer 
to perform the trade-offs among them. 

CONDUIT is built on top of the MATLAB/SIMULINK system modeling and 
analysis environment and the C-0 computer-aided, optimization-based para- 
metric design software package. It includes a graphical block diagram editor, 
a graphical spec editor, a spec library, an editor for the initialization file and 
multiple layers of supporting analysis plots. This avoids the manipulation of 
MATLAB “m” files or C-0 C-type files. The user only has to select and drag 
and the system automatically updates the relevant scheme, spec or plot. 

CONDUIT has two modes of operation: setup and run. 


1.1 Setup Mode in CONDUIT 

The setup of a problem in CONDUIT consists of three parts: 
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1. First the designer has to define (or import) into SIMULINK a simulation 
of the aircraft and the controller. The real aircraft control system consists 
of many controller loops, all of them controlling the high-order, nonlinear 
aircraft dynamics. Usually people do not use the nonlinear dynamics in 
design because it is too complex and needs special software for simula- 
tion. Instead they use linearized dynamics. CONDUIT uses SIMULINK 
for simulation, so that the block diagram can include nonlinearities in the 
dynamics. The fact that some nonlinearities, especially saturations and 
dead-zones, are important to aircraft control system design is one moti- 
vation for the CONDUIT design approach. Designs based on a purely 
linear model of the controlled aircraft that ignore saturation, delays, and 
dead-zone may be overly optimistic. 

There usually exist several models for the same system, corresponding 
to particular flight conditions (speed: hover, low speed, forward flight or 
mission tasks: precision tasks, aggressive task, etc.). For a thorough ex- 
amination of the problem, the control systems of all these models should 
be investigated and optimized. 

The control law model must include the design parameters such as gains 
and time constants, as well as the inputs and outputs of the system. The 
inputs are used to apply test input signals. The outputs are used to eval- 
uate system responses. 

2. The next step for the designer is to define the set of design specifications. 
To help the designer, there are five graphical libraries comprising over 50 
specifications in CONDUIT. These are organized into the following sub- 
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libraries: 


• Generic Specifications 

• Fixed-Wing Longitudinal 

• Fixed-Wing Lateral/Directional 

• Rotor Hover Low Speed 

• STOVL (Short Take-Off / Vertical Landing) 

Most of the specifications belong to some aeronautical design standards: 
rotorcraft specifications belong to ADS-33D [19] and aircraft specifications 
belong to MIL-STD-1797 [20] and Mil STD 9490. 

The designer only has to select the appropriate design specification from 
the available libraries, drag them into the HQ Window (Handling Qualities 
Window), and configure the simulation appropriately for each specification 
using the HQ Editor (Handling Qualities Editor). 

All the specs are created using the same color map: the Level 1 region is 
blue; the Level 2 region is magenta; the Level 3 region is red. The splines 
that define the boundaries between levels can be graphically altered to 
design new specifications. 

To “wire” a specification to the SIMULINK simulation model the designer 
has to use the graphical “HQ editor” . Here, the user declares each speci- 
fication to belong to one of the following five classes: 1) hard constraint, 2) 
soft constraint, 3) objective, 4) summed objective, or 5) check only. Those 
specifications that are declared to be “summed objective” are treated as 
soft constraints until the design reaches Level 1. Then, all the summed 
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objective specs are added together and treated as a single objective. At 
this stage, an improvement in any component of the summed objective 
improves the design s performance. The selection of the specification class 
defines the problem for the C-0 optimization process. In aircraft design, 
examples of choosing the specifications are the following: stability specifi- 
cations are hard constraints; most of the other specifications are soft con- 
straints (bandwidth and time delay (frequency-domain point spec), quick- 
ness (time-domain point spec), attitude hold (time-domain envelope spec), 
interaxis coupling (time-domain point spec and envelope spec)); crossover 
frequency, actuator energy, and actuator saturation are usually perfor- 
mance objectives. 

The input and output port connections for each specification are indicated 
in an information box in the spec editor. 

In case the designer needs to add a specification that is missing from the 
library he can create it and add it to the library with the “spec maker” 
facility. 

3. The designer must also set up a small initialization file to define problem- 
dependent constants such as the simulation time-step, the performance 
precision, denoted by epsilonF, and the test input signals. 

In the following paragraphs some details regarding the problem setup will 
be presented. “Stability specification” and “Gain margin and phase margin 
specifications” refer to the specs that must be included in the HQ Window in 
order to insure a stable system. “Simulation time-step” and “Gradient step-size” 
refer to some options that need to be specified in the initialization file. 
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1.1.1 Stability Specifications 

Stability criteria are the most important specifications for aircraft control sys- 
tem design and therefore are the very first selected in the performance chart. 
The “Generic Specifications” section of the CONDUIT spec library contains two 
specs designated to evaluate the system stability: EigLcGl (eigenvalue test) and 
StbMgGl(gain and phase margin test). Usually these two specs should be used 
together. Using only the eigenvalue test, it is not possible to determine the clas- 
sical stability margins that define robust stability. However, the eigenvalues are 
a definitive test of closed-loop stability. In contrast, the classical gain and phase 
margins are insufficient to determine stability. In part this is because they apply 
only to specific loops, rather than the complete closed-loop system. This is also 
due to errors in MATLAB. MATLAB computes gain and phase margin based 
on a Bode diagram without first checking the stability. This makes it possible 
that unstable systems appear to have good gain and phase margin, as in the 
following example. 

Example 1: Take the system in Figure 1.1. Checking in MATLAB: 

>> [g,p] = bode(num, den, logspace(— 1, 5, 1000)) 

>> margin(g, p, logspace(— 1, 5, 1000)) 

The MATLAB function “margin” shows gm=14.18 and pm=45.05, as shown in 
Figure 1.2, which are good stability margins. The poles of the closed-loop system 
are: 

-10.5415, -0. 7551+2. 2655i, -0.7551-2.2655i and 0.0516, so the system is unstable. 
The Nyquist plot constructed for this system, with the contour modified to take 
a small detour around the pole, at s=0, is shown in Figure 1.3. It is apparent 
that the polar plot of the open-loop system encircles the -1+jO point once as u> 
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Transfer Fcnl 


Figure 1.1: Example of unstable closed-loop system for which MATLAB indi- 
cates good stability margins 




Figure 1.2: Gain margin and phase margin for the system in Fig. 1.1 evaluated 
by the MATLAB function “margin” 
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Figure 1.3: Nyquist plot for the system of Fig. 1.1 


varies from — oo to +oo and thus this closed-loop system is unstable. 

The solution to eliminate the instability is to use the eigenvalue test. It is 
recommended to use EigLcGl and StbMgGl both as hard constraints, to be 
evaluated and optimized in Phase 1. 

1.1.2 Gain Margin and Phase Margin Specification 

There are some stable systems that have both a gain margin increase and a gain 
reduction margin, i.e. both the gain increase and gain decrease must be limited. 
CONDUIT allows one value for the gain margin, so in this case it will take the 
smallest absolute value of the two of them as the gain margin. This means that 
the system can increase or decrease by that gain margin and still remain stable. 
Example 2: 

The antenna problem, part III, chapter 3, has two gain margins for the following 
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parameters, called the nominal point: 


dpp.l = 1.08000e + 02 dppj = 5.00000e + 02 

dppJ. = 5.40300e + 01 dpp. 8 = 5.00132e + 03 

dpp. 3 = 2.62000e + 03 dpp. 9 = 7.54000e + 02 

dpp A = 1.57100e + 03 dpp. 10 = 2.00000e - 01 

dpp. 5 = 1.41370e + 00 dpp. 11 = 6.67500e - 01 

dpp. 6 = 4.00000e + 00 

This system is stable, having all poles in the left half plane. Figure 1.4, Figure 
1.5, and Figure 1.6 show that the system is conditionally stable. The polar plot 
of the open-loop transfer function cuts the real axis at -0.4636 and at -48.4172. 
This shows that the 

gain can increase by 1/0.4638 (20*logl0(l/0.4638)dB = 6.67dB) and the 
gain can reduce by 48.4172 (20*logl0(l/48.4172)dB = -33.7dB) 
before instability occurs and defines both a gain increase margin and a gain 
reduction margin. Figure 1.7 shows the gain reduction margin. 

For this case CONDUIT will take max(abs(6.67),abs(-33.7))dB=6.67dB as 
the gain margin. CONDUIT uses its own “margin” function that fixes the prob- 
lems with MATLAB’s “margin” function. MATLAB’s “margin” function doesn’t 
show all the gain margins and phase margins. 

1.1.3 Simulation Time-Step 

The real aircraft control system consists of many digital controller loops, with 
the inner controller loops faster than the outer controller loops, all of them 
controlling the high order, nonlinear aircraft dynamics. Designers prefer to work 
with continuous systems because they are easily understandable. Simulations, 
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Figure 1.4: Nyquist plot for Example 2 as obtained from MATLAB without 
specifying the frequency range 



Figure 1.5: Nyquist plot for Example 2 on the [10 11 , 10 3 ] frequency interval 
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Figure 1.6: Nyquist plot for Example 2 on the [10 2 , 10 4 ] frequency interval 



Figure 1.7: Bode plot for Example 2 
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such as those produced by SIMULINK, run faster using discrete systems. Then 
there is a choice to use in simulation: 

1. Use a continuous model and continuous controller 

This has the advantage of fidelity to the real aircraft - aircraft fly in con- 
tinuous time. Also, most designers have more experience with continuous- 
time systems and better intuition. Lastly, the specifications are based on 
continuous-time behavior. 

However, the simulation of continuous-time systems is relatively slow. 
Moreover, the actual implementation of the controllers is in discrete time. 

2. Use a continuous model and digital controllers, each with its own sample 
rate 

This has the advantage of very accurate results because of the use of real- 
istic digital controllers, but slow simulation. 

3. Use a continuous model and continuous controllers, discretized at a single 
sample rate 

This has the advantage of choice 1 as well as a relatively fast simulation. 
The disadvantage is that the solution obtained is less accurate. However, 
it is easy to obtain sufficient accuracy by appropriate choice of the step 
size. 

CONDUIT uses the third option. The designer has to provide a continuous 
model and continuous controllers for the block diagram. The designer must also 
choose a sampling interval T, for the simulation. The value of T is placed in the 
initialization file. Some general rules of choosing T are: 
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1. An absolute lower bound for the sample rate (1/T) is stated by the Shannon 
sampling theorem [12]. If it is desired that the closed-loop system have a 
certain bandwidth, say u b , then the sample rate must be at least twice this 
required closed-loop bandwidth. 


2. If white noise disturbances are the dominant source of error in the system, 
then the sample rate should be faster. Theoretically, all real signals have 
spectral content at all frequencies so there is always information lost due 
to sampling. By [12], testing the RMS (root mean squared) value of the 
steady-state value of the covariance of one output, y ^Efy 2 ], for the same 
design and design parameters, but different sample rates, shows that for 
good random disturbance attenuation a sampling interval of T = -2- 

° 10o; 6 

would be a good choice. The relative errors grow quickly when sampling 
slower than this multiple, whereas very little can be gained by sampling 
faster. Thus, 


7 r _ 7T 

< T < — 

lOoij, u b 


( 1 . 1 ) 


CONDUIT proposes that the sample rate ± be 20 times bigger than the 
fastest mode of the closed-loop system, determined from the eigenvalues of the 
closed-loop system. Usually the fastest mode for aircraft lies between 2Hz and 
10Hz, so the sampling interval, T, is between .005 secs and .25 secs. 

The optimal controllers computed by CONDUIT are based on a discretized 
version of the aircraft and its controller. The discretization uses sampling interval 
T. the controller actually used by CONDUIT is shown in the dtest.m block 
diagram in CONDUIT. 

The designer should be aware that for different values of T in the initialization 
file, CONDUIT can provide slightly different controllers. 
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1.1.4 Gradient Step-Size 

C-0 evaluates the gradients of the optimization specs at every iteration during 
the optimization process. The increment used to compute gradients in theoretical 
problems for which the data is perfectly accurate, is 

ppdelta = maxirteps * max( I.eO, fabs(x[i])), pdeltal) (1.2) 

where x[z] is the parameter perturbed in computing the gradient, fabs is the 
absolute value, rteps is the square root of machine precision, pdeltal is an arbi- 
trary constant step, usually 0. This formula shows that the gradient step-size is 
usually rteps * fabs(x[i}). It is rteps if x\i) is smaller than 1 and it is pdeltal if 
the user chooses a bigger step. 

For real problems, such as those optimized by CONDUIT, the increment used 
is 


\J epsilonF * max(l.e0, /ahs(x[i])) (1.3) 

where epsilonF is chosen by the designer and represents the precision of the 
optimization specs (C-0 allows different epsilonF between specs). This way the 
gradient step-size is usually \J epsUonF * fabs(x[i]) and it is y/epsilonF if x[z] is 
smaller than 1. The advantage of using \/ epsilonF *fabs(x[i ]) instead of pdeltal 
is that this bigger step is proportional to the parameter, not a constant value. 
Moreover it can be set up from the initialization file. The current implementation 
of CONDUIT, to simplify the designer’s job, allows only one epsilonF for all 
the specs. The designer has to set the value in the initialization file. The default 
value is 0.0002 and was chosen following a set of tests on aircraft optimization 
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problems [24]. Note that this is substantially larger than the increment used in 
theoretical problem. 

1.1.5 Creating new specs with the “Specmaker” 

The Specmaker facility gives flexibility to the designer to create and use any 
performance criterion he wants in CONDUIT. There are 5 types of specs that 
can be added to CONDUIT [25]. 

• time point 

• frequency point 

• time line 

• frequency line 

• LOES (lower-order equivalent system) 

For time point and frequency point specs the performance criterion must 
consist of 2 lines that separate the 3 levels of performance. For time line and 
frequency line specs the performance criterion must consist of 2 envelopes ( closed 
curves) that encircle the Level 1 and Level 2 regions, respectively. Outside both 
envelopes is Level 3. Examples are given in the following chapters. 

Specs based on lower order equivalent systems (LOES specs) are common in 
aricraft control. The basic idea is as follows. Given a time or frequency response 
plot obtained either experimentally or from a detailed high-order simulation. 
Approximate this given response by the response of a lower-order and much 
simpler system. The lower order system is characterized by a relatively small 
number of parameters. These parameters are chosen to minimize some precise 
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measure of the difference between the lower order model and the real system. 
The actual specification is given in terms of the lower order model. 

1.2 Run Mode in CONDUIT 

Once the problem is set up the second phase, design tuning, follows. During 
this phase the designer cannot change the simulation model, performance speci- 
fications, or initialization file. But he can change the current design parameters, 
evaluate the specifications for every set of design parameters, and use the analy- 
sis functions of the system. And CONDUIT will run C-0 to optimize the design 
parameters. 

A distance algorithm in CONDUIT translates the location of the design point 
(in case there is a simple specification) or the location of the worst design point 
(in case there is a functional specification) on each of the graphical specification 
criteria to a numerical rating. All of C-O’s “good-bad” factors are scaled - the 
scaled good values become 1 and the scaled bad values become 2. A rating of 
“1” for the scaled design point indicates that it lies on the Level 1/2 border. A 
rating of “2” indicates that the design point lies on the Level 2/3 border. These 
numerical ratings are used by C-0 to tune the design and they appear in the 
Pcomb (Constraint Table). CONDUIT feeds the specifications in the form of 
the “problem description file” to C-O. 

C-0 implements multicriterion parametric optimization. It uses FSQP (Fea- 
sible Sequential Quadratic Programming) [1] to solve the following general opti- 
mization problem: 

minobjf(x) Vi (1.4) 
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subject to : soft?(x) < 0 Vj 

hardj(x) <0 \/k 

hard-boundi(x) < 0 V/ 

where objf, soft®, hardjj. are the scaled values of objectives, soft constraints, hard 
constraints (except hard bounds), respectively and where hard.boundi represents 
hard bounds, an important special case of hard constraints. 

• An objective is a specification of a quantity that should be optimized (mini- 
mized or maximized). Typically, multiple competing objectives are present. 

• A hard constraint is a specification of a quantity that must achieve a spec- 
ified threshold. A design for which a hard constraint is not satisfied cannot 
be satisfactory. 

• A soft constraint is a specification of a quantity that should achieve , or at 
least approach, a specified threshold, i.e., should be optimized as long as 
this threshold is not achieved. 

• A hard bound is a constraint of the form a < dpi < b, where a, b are 
constants and dpi is the i th design parameter. 

C-0 uses the min/max optimization criterion: 


nun (maxi<i< m a:jfj(x)), 0 < real (1.5) 

where fj(x) is the ith specification and the ot\ are user-specified weighting coeffi- 
cients. The advantage of this formulation is that the optimal value of x can be 
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placed anywhere in the region of the parameter space bounded by the minima 
of the individual criteria by appropriate choice of the a*. 

The optimization process is divided into three different phases. 

Phase 1: 

Not all hard constraints are satisfied. Equation 1.4 takes the form 


min maxi<hard£ (x) 

(1.6) 

subject to : hard .bound) (x) <0 VI 

(1.7) 


Phase 2: 


All hard constraints are satisfied (have values < 1). Not all objectives and 
soft constraints are Level 1. Equation 1.4 takes the form 


min maxj j{objf (x), soft® (x)} 


(1.8) 

subject to : hard^(x) < 0 

Vfc 

(1.9) 

hard.boundi(x) < 0 

VZ 

(1.10) 

Phase 3: 

All hard constraints and soft constraints are satisfied and all objectives are 

Level 1. Equation 1.4 takes the form 

min maxjobjf(x) 


(1.11) 

subject to : soft®(x) < 0 

Vj 

(1.12) 

hard^(x) < 0 

VJfc 

(1.13) 

hard-bound) (x) < 0 

VZ 

(1.14) 
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The active objective (constraint), i.e. the one C-0 is working on, is displayed 
in the title of the Pcomb chart in C-0 and has a name according to the opti- 
mization phase. In Phase 1 it is named MAX_HARD, in Phase 2 it is named 
MAX_COST_SOFT, and in phase 3 it is named MAX-COST. The value of the 
displayed active objective (constraint) is scaled so that 1 shows the Level 1 /2 bor- 
der and 2 shows the Level 2/3 border. MAX-HARD and MAX_COST_SOFT 
are always positive, because the active constraint is Level 3 or Level 2, while 
MAX-COST is always negative because the active objective is Level 1. 
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Part III 


Specific Work 
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Chapter 2 


F-14 Problem 

2.1 Stability and Control Augmentation Sys- 
tems 

Most modern commercial and military aircraft have Stability and Control Aug- 
mentation Systems (SCAS) in order to meet the flying qualities requirements. 
SCAS decrease the pilot workload by controlling the position and angular rate 
of the aircraft to match the values of pilot input references. The system also 
reduces the coupling that normally exists between different pilot inputs. The 
SCAS outputs add to the pilot’s. SCAS outputs move only the control surfaces 
and not the cockpit controls. Automatic pilots, which replace the human pilot 
when they are in use, are expected to move the cockpit controls. 

In the case of high-performance military aircraft, where the pilot may have to 
maneuver the aircraft to its performance limits and simultaneously perform tasks 
such as precision tracking of targets, specialized control augmentation systems 
are needed in addition to the stability augmentation systems [ 13 ]. 

A pitch-axis control augmentation system is a specialized SCAS and it is 


27 



used in several inodes of operation corresponding to the aircraft type and the 
required performance. 

For a fighter aircraft, where high maneuverability or “agility” counts the 
most, the suitable controlled variable for the pitch axis is normal acceleration 
(or load factor) on the aircraft. This is the component of acceleration measured 
by an accelerometer in the negative direction of the aircraft z-axis [14]. It is 
directly relevant to performing a maximum-rate turn and must be controllable 
up to the structural limits of the airframe, or the pilot’s physical limits. 

A second common mode of operation for a pitch-axis control augmentation 
system is as a pitch-rate command system. Control of pitch rate is used when a 
situation requires precision tracking of a target and is also preferred for approach 
and landing. 

These two modes are contradictory. A control system that has a good normal 
acceleration step response may have a pitch-rate response with a very large 
overshoot, and conversely, a reduction in the pitch-rate overshoot may lead to a 
sluggish normal acceleration response. For fighter aircraft air-combat modes a 
pitch-rate overshoot is required for good gross acquisition of targets and a dead- 
beat pitch-rate response is required for good fine tracking. Sometimes the two 
control schemes are blended together to give the pilot control over pitch rate at 
low speed and normal acceleration at high speed. 

There exists also a third mode of operation, and this is as an attitude (or 
angle-of-attack) command system. This is the one used for the Grumman F-14 
aircraft and will be presented in the following sections. 
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2.2 Description 


In one functional mode, the F-14 uses an angle-of-attack command system for 
the pitch-axis control augmentation system. 

The F-14 problem is to design a system which controls the pitch response 
to pilot stick commands. The controller performance is tested for speed, angle- 
of-attack tracking when following an a-command, and for disturbance rejection. 
The stability, the magnitude of the tail rate, the derivative of the normal accel- 
eration, and the shape of the pitch rate are also controlled. 

The problem is derived from a master’s thesis example [5]. Some transfor- 
mations changed the block diagram into an equivalent one, Figure 2.1, that 
matches the Grumman F-14 Benchmark Control Problem formulations [8], [9], 
[10]. These transformations are: 

The pitch-rate lead filter was changed from to s+d P - 2 , 

° s— dp - 1 s+dp-1 

The stick prefilter was changed from — j-1 — - to -r-4 — r. 

3^5+1 dp-73+1 

The a-sensor low- pass filter was filter changed from — r-1 — to -r-k — - 

5F3 ,+ 1 <*P-3s+l 

where 

dpi = negative of the pitch-rate lead filter pole location. 

dp 2 = negative of the pitch- rate lead filter zero location. 

dp 3 = negative of the inverse of the a-sensor low- pass filter pole location. 

dp 7 = negative of the inverse of the stick prefilter pole location. 

The controller gains were changed from 



(2.1) 

( 2 . 2 ) 
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This implies that dpJa , appears in the forward path multiplying the output of 
the stick prefilter but has no effect on the pitch-rate and a-feedback gains where 
dp 4 = pitch-rate feedback gain 
dp 5 = a-feedback gain 
dp$ = gain of the filtered command signal 
dps = common integral gain 

White noise filtered through the Dryden wind gust model was added. For 
optimization the EigLcGl, RMSlspc and BnwPiL2 specs were added. The spec- 
ifications are shown in Figure 2.6. 

The system has to have the form of Figure 2.1, where parameters dp A through 
dp . 8 are to be selected to meet the performance objectives and constraints listed 
in section 2.4. Selecting the gains dpA through dp . 8 to yield good closed-loop 
response to a step input at input 1 and white noise at input 2 corresponds to a 
multi-input multi-output design problem. 

2.3 Block Diagram 

Figure 2.1, Figure 2.2, and Figure 2.3 show the block diagram of a pitch-axis con- 
trol augmentation system. The aircraft dynamics represent those of the Grum- 
man Aerospace F-14 flying at 35,000 ft in level flight at a velocity of 690 ft/sec. 
The model is linearized, time invariant and continuous time. 

The components of the model are: 

• short-period longitudinal aircraft dynamics 

• tail surface actuator 
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Figure 2.1: Grumman F-14 pitch-axis control system 
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Figure 2.2: Alpha response model 
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• control system, consisting of a stick prefilter, an error integrator, pitch-rate 
lead filter, and an angle-of-attack filter 

• Dryden wind-gust model 


An error integrator has been included to make the control system Type 1 , thus 
ensuring that the aircraft will hold a zero angle-of-attack trajectory when no 
pressure is applied to the control stick. The Dryden Model filters white noise to 
form the vertical gust velocity, u gu „fi which is then filtered further to give the 
angular rate gust component, q gus t- 

The fixed parameter values shown in Figure 2.1 are: 

Tail Servo 



c = 0.707 
o> 0 = 2.49^ 

u sec 

Grumman also supplied values for the design parameters dp.l through dp. 8. 
These were the results of their quadratic-based optimization system, CASCADE. 
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These values of: 


dp A = 4.14400e + OOg dp. 5 = 6.77000e - 01 

dpJ2 = 2.97100e 4- 00 g dp Jo = -1.74600e + 00 

dp- 3 = 3.95900e - Olsec dp. 7 = l.OOOOOe - Olsec 

dp A = 8.15600e - 01 dp.8 = -3.86400e + 00 

were used as the initial values for the independent design parameters dpA 
through dp. 8. 

2.4 Performance Objectives and Constraints 

The second major component, after the system to be controlled and the con- 
troller structure, of a control system design problem is the set of specifications. 
CONDUIT provides a very large collection of aircraft control system specifica- 
tions in its on-line specs catalog. All the specifications are illustrated in Fig. 
2 . 6 . 

Specifications are described in an abbreviated fashion. The purpose of the 
spec is given general terms. This is followed by the CONDUIT description of 
the spec, written in bold text and followed by the word, “code.” The following 
lines are written in pseudo-code, with comments. The first line after the header 
indicates whether the spec applies to the open-loop or closed-loop system. This 
same line specifies the inputs and outputs used in computing the spec. The 
next line, in the absence of comments, designates the spec as on objective, hard 
constraint or soft constraint. This line also indicates the range of frequencies or 
times over which the spec is evaluated. The last line gives the goal the spec is 
intended to achieve. Pictures of all the specs can be found in Fig. 2.6. 
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Performance Objective 

This specification, when used as a performance objective, generally tries to 
make the closed-loop bandwidth larger and the time delay smaller. 

1 . The pitch bandwidth specification requires a lower limit for pitch attitude 
bandwidth and an upper limit for time delay. Aircraft maneuvers are 
classified in different flight phase categories, corresponding to the required 
speed of the maneuvers, precision of the tracking, precision of the flight- 
path control [20]. This spec applies to all categories so it is specified as 
Categories A&D. 

BnwPiL2 code: 

closed-loop, input 1 (the pilot pitch command), output 7 (a (pitch re- 
sponse)). 

objective, u G [0.1,20] Hz 

maximize bandwidth and minimize time delay 

Constraints 

1. The first stability criterion of the total feedback system (i.e., with both 
feedback paths active) requires nonpositive eigenvalues: 

EigLcGl code: 

closed-loop, from any input to any output 

[a,b,c,d]=closed loop model, unstab equals sum of positive eig(a) if there 
are some positive eigenvalues or unstab equals max of negative eig(a) if all 
eigenvalues are negative 
hard constraint 

minimize unstab : Level 1/2 = 0.001, Level 2/3 = 0.002 
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2. The second stability criterion of the system is: 
gain margin > 6 [dB] 

phase margin > 45 [deg] 

StbMgGl code : 

open-loop, input 3 broken-loop input, output 6 broken-loop output 
hard constraint, u e [0.1, 100] Hz 
maximize gain : Level 1/2 = 6 [dB], Level 2/3 = 4 [dB] 
maximize phase: Level 1/2 = 45 [deg], Level 2/3 = 35 [deg] 

3. The angle of attack, a [deg], response to a pilot-step input of 2.0 [deg] is 
required to match that of a given critically damped second order system. 
The resulting difference is to be minimized with an acceptable value being 
0.2 [deg]. Because the system is linear this is equivalent to acceptable 
design for a error less than 10%, satisfactory design for a error less than 
5%. 

aerrspc code: 

closed-loop, input 1 the pilot pitch command, output 2, q=pitch rate 
soft constraint, t e [0.1,3] sec 

minimize a error, Level 1/2 = [-0.05 0.05], Level 2/3 = [-0.1 0.1] 

4. The maximum velocity of the tail surface deflection servo must be less 
than 25.0 [deg/sec] for the 2.0 [deg] pilot-step input. Because the system 
is linear, this is equivalent to satisfactory design for ratio between tail rate 
and command less than 12.5%. 

tailspc code: 

closed-loop, input 1 the pilot pitch command, output 1 the tail rate 
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soft constraint, t G [0.05, 3] sec 

minimize tail rate, Level 1/2 = [-12.5 12.5], Level 2/3 = [-15 15] 

5. The rate of change of acceleration (i.e. jerk) experienced by the pilot in 
the vertical direction must be greater than zero during the response to the 
pilot-step input. The normal acceleration is defined as 
normal acceleration = Uo*q + L*q — u, where Uq appears in Figure 2.1, 
at aircraft dynamics, q is the pitch rate, q is the derivative of the pitch 
rate, u> is the output from Sum6 in Figure 2.1, and L = 22.8 [ft]. The 
derivative of normal acceleration is obtained from the normal acceleration 
by a multiplication by in the block diagram. The derivative of normal 
acceleration is scaled by the acceleration due to gravity, g 0 = 32.2^. 

For this spec two signals were considered, Signal 1 corresponding to an 
input of 2 [deg] and Signal 2 corresponding to an input of 5.72 [deg]. 

This constraint must be positive, but the value does not matter. Thus 
a saturation was introduced in the block diagram, with limits ±1. If the 
value is smaller than -1, the spec will show -1, which is Level 3. If the value 
is bigger than 1 , the spec will show 1 and is Level 1 . If the saturation were 
omitted then the envelope spec would need an expanded range for Level 
1, and this would decrease the readability of the spec, 
jerkspc code: 

closed-loop, input 1 pilot pitch command, output 4 the jerk ,a = 
acceleration) experienced by the pilot 
soft constraint, t G [0.1, 1.5] sec 

maximize jerk, Level 1/2 = [0 1] [1/sec], Level 2/3 = [-0.2 1.1] [1/sec], 
saturation [-1, 1] 
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6. The pitch rate (q) [rad/sec] response to stick inputs [rad] is required to fit 
the lower-order system [6]: 


ir(1.7s-H)e~ TJ 

S 2 + 2^U)qS + U)q 


within the bounds of 
2.122 < cj 0 < 3.5214 
0.5 < C < 0.707 
r < 0.1 

for 0.2 < u < 20 

Because K is unconstrained, this requirement constrains only the transient 
response. Equation 2.3 and the associated bounds can be translated into 
an equivalent gain and phase requirement. 

The gain requirement is: 

y/T^T{L7uj)i g a i n ^1 + (1.7a;) 2 

min — •- = < < m ax 

s> 2 dC9am ~ 

(2.4) 

for 0.2 < u> < 20 

The minimum is attained for u 0 = 2.122 and C = 0.707. The maximum is 
attained for uo = 3.5214 and £ = 0.5. Equation 2.4 is equivalent to 

v/ 1 + (»-*■>>’ < gain , 

v'(iMife)TT(270707^ “ *«“*" - 

(2.5) 

jA + ( 1.7a>) 2 

" (l*4) ! ) 2 + (2*0'55sn) 2 

for 0.2 < u < 20 
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The phase requirement (time delay) is: 


(j) 

2^ 360 

min (arctan 1.7 u — arctan ^ ra;)— — [deg] < < 

J _ ^ Z7T 

CJo 


(2.6) 


U) 

, 2< ^ 360, , , 

max (arctan 1.7u — arctan A to;) — — [deg] 

^o.C.t 1 _ ( ) 2 27r 

U^o 


for 0.2 < u> < 20 


The minimum is attained for Uq = 2.122, £ = 0.707, and r = 0.1. The 
maximum is attained for u>o = 3.5214, ( = 0.5, and r = 0. Equation 2.6 is 
equivalent to 


u> 


(arctan 1.7a; — arctan 


2 * 0.707 ocn 

- 01 * u; )^r[ de g] ^ ^ 

1 - ( ) 2 27r 

v 2.122 ; 


2*0.5„ -~z - 360 r 


(2.7) 


(arctan 1.7u; — arctan ^-5214 )— — [deg] 

1 - ( ) 2 2?r 

v 3.5214' 


for 0.2 < to < 20 

Another form of this requirement is the “YF-17” frequency-response crite- 
rion [16]. It shows the matching of the high order system (HOS) pitch rate 
to stick input with some gain and phase curves. This criterion includes a 
constraint for the degain. Figures 2.4 and 2.5 show a comparison between 
( dcgain*LOS ) and “YF-17” criterion. For our scheme dcgain=0.5697 at 
the nominal point. Therefore gain and phase curves of 

0.5697 * ^\ 7s — and “YF-17”. will be compared. 

s 2C s 

-5 + — + 1 

U)q u ) 0 
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Gain [dB] 



Frequency [rad/sec] 


Figure 2.4: Comparison between YF-17 criterion and analytic curve of dc- 
gain*LOS. This is the gain spec. The YF-17 Level 1 region is indicated in black. 
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YF-1 7 Phase Criterion (pitch axis) 



Frequency [rad/sec] 


Figure 2.5: Comparison between YF-1 7 criterion and analytic curve of dc- 
gain*LOS. This is the phase spec. The YF-17 Level 1 region is indicated in 
black. 
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The theoretical criterion, which considers the gain and phase limits for the 
pitch rate (divided by degain), will be used in the optimization. The prac- 
tical criterion YF-17” will be shown as “check only” in the performance 
chart, without any role in optimization. The degain option for gain is set 
to 20*logl0(0.5697)=-4.8871 [dB], 

YFGAspc and YFPHspc code: 

closed-loop, input 1 pilot pitch command, output 3 pitch rate 
soft constraint, u € [0.2, 20] Hz 

7. RMS g’s at the pilot caused by turbulence have to be minimized. Distur- 
bance rejection is an important aspect of any control system, so there is 
a soft constraint which takes care of errors due to random disturbances. 
Berman and Gran [13] suggest that the sample rate selection for aircraft 
autopilots should be done based on its effect on disturbance rejection and 
discuss an application of digital design and sample rate determination to 
pitch-plane control of the Grumman V/STOL Design 607A. 

Disturbances enter a system with various characteristics ranging from steps 
to white noise. For purposes of determination of sample rate, the higher 
frequency random disturbances are the most influential; therefore, we will 
concentrate on their effect. In other words, disturbances that are fast 
compared to the plant and the sample rate will be compared. That is, 
where the noise can be considered to be white. 

In case the control system uses a good continuous controller, the magnitude 
of the error response represents a lower bound on the magnitude of the 
error response that can be hoped for when implementing the controller 
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digitally. Some degradation from the continuous design occurs because of 
the sampling. In order to analyze the degradation of the digital controller 
as compared to the continuous controller, it is important to consider the 
effect of the noise consistently, with both the continuous and the digital 
controllers. 

Suppose the system is continuous and represented by 

x = Ax + Bw (2.8) 

where system matrices A and B describe a closed-loop system including a 
continuous controller and w = [wiW 2 -..w n ] T is the noise. Denote the power 
spectral density of w as Rwpsd (alternatively referred to as the “white-noise 
intensity” or “mean-square spectral density”). Then the covariance matrix 
of w is E[w(t)w T (t + r)] = Rwpsa^r). 

This assumes that the disturbance w(t) is white noise, which means that 
w(t) and w(s) are uncorrelated for all t ^ s. 

Typically, if there is more than one process noise component (i > 1), 
one has no information on the cross-correlation of the noise elements and 
therefore R wp sd is selected as a diagonal matrix. 

The steady state value of the covariance matrix of x is given by the Lya- 
punov equation 

AX + XA t + BR wp8d B T = 0 (2.9) 

The solution to this equation, X=E[x(t)x T (t)]) represents the amplitude 
of the random response of the state due to the excitation from w. X 
can be used to establish a baseline against which controllers are com- 
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pared. For this particular problem the root mean square (RMS) value, 
computed as the square root of the steady state value of the covariance of 
y> \j (E[Cx(t)x T (t)C T ]) = yj (CXC T ) will be used, y is a scalar because one 
particular output is considered, e.g. pilot g’s. And R wpsd is assumed to be 
1. 

The maximum RMS(root mean squared) must be less than 0.04 [8,9], 

RMSlspc code: 

closed-loop, input 2 white noise input to the Dryden wind gust model, 
output 5 pitch acceleration (g’s) 
soft constraint 

minimize RMS, Level 1/2 = 0.04, Level 2/3 = 0.1 where 
X = lyap(A,B(:,2)*B(:,2) T ) 
gVar = C(5,:)*X*C(5,:) T 
RMS = yj (gVar) 

It is also necessary to evaluate X when the system has a digital controller 
for the identical excitation applied to the plant. Suppose the discrete 
equivalent system is 

x(k + 1) = A d x(k) -f B d w(k) (2.10) 

where the process noise w(k) is a random sequence with zero mean, that 
is, E[w(k)] = 0, and has no time correlation or is “white noise”, 
E[w(i)w T (j)] = 0 if i ^ j, and has covariances or “noise levels” defined by 
E[w(k)w T (k)] = R w . The desired result, called the discrete Lyapunov equa- 
tion, is 


A d XAj + G = X 


( 2 . 11 ) 



where G = / 0 T A d (r)BR wpsd B T Ad (r)dr. If T, the sampling interval, is 
much shorter than all system time constants, that is, Ad = I and Bd = BT, 
then the integral is approximately G = Bd R ^ p<j Bj = B d RwB d [12]. In or- 
der to evaluate the effect of sample rate on the performance of a controller 
in the presence of white plant disturbances, first Eq.2.9 should be evaluated 
to find the baseline covariance X and then Eq.2.11 should be repeatedly 
evaluated with varying sample rates to establish the degradation versus 
sampling. The RMS value is the quantity that is typically measured. In 
the discrete case, for the F-14 particular problem, the RMS value is com- 
puted in MATLAB as follows: 

X = dlyap(A d , B d (:, 2) * B d (:, 2 ) t /T) 
gVar = C d (5, :) * X * C d (5, :) T 
RMS = yj (gVar) 

Taking relatively fast sampling (twenty times the fastest frequency mode), 
the RMS value is almost the same for both the continuous and digitized 
system, so that either one can be used, in RMS computation. In case that 
dt=l/50 and nominal design parameters are used rms con tinuous = 0.0367244 
and 

rms discre te = 0.0367230. 

2.5 Setup of the Initialization File 

Looking at the closed-loop eigenvalues at the nominal starting point, the fastest 
pair of eigenvalues is -9.84 ± 9.57j, giving the largest u = 13.72 and the 
largest v = 2.18 Hz. The sampling rate is chosen to be 50 samples per second. 
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epsilonF was set at 0.0002 (the default value) and 0 (machine precision), and 
tests were made for both values. 

The input signals InpSigl and InpSig2 are step signals, corresponding to 2 
and 5.72 [degjrespectively. 


2.6 Optimization 

The F-14 optimization starting at two different initial points will be presented: 
a satisfactory starting point found by Grumman using their quadratic-based 
optimization system CASCADE, called the nominal point, and an arbitrarily 
chosen point, called the trivial point given below: 

dp A = 1 dp . 5 = 1 

dp_2 = 1 dp . 6 = —1 

dp. 3 = 1 dp. 7 = 1 

dp. 4 = 1 dp. 8 = -1 

In both cases CONDUIT was able to improve the performance and to deliver a 
satisfactory result. The optimization results using two values for epsilonF are 
shown in Table 2.1. 

The optimal point achieved when the optimization started in the nominal 
point and used epsilonF = 0 was: 

dp. 1 = 4.14184e + 00 dp. 5 = 6.77324e - 01 

dp.2 = 2.97359e + 00 dp. 6 = -1.74871e + 00 

dp. 3 = 3.93938e - 01 dp. 7 = 5.28613e - 02 

dp_4 = 8.22272e - 01 dp. 8 = -3.86422e + 00 
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Table 2.1: Comparison between the optimal points obtained by CONDUIT for 

the F-14 problem == =^==^=^===^=:^=== 

Starting point Optimal using Optimal using 

epsilonF=0 epsilonF=0.0002 

Nominal MAX. COST M AX .COST M AX .COST 

-0.503495 -1.2329 -0.955191 

Trivial MAX. COST JSOFT MAX. COST JSOFT MAX.COST^OFT 

19.6509 0.0779831 0.376201 

The optimal point achieved when the optimization started in the nominal 
point and used epsilonF = 0.0002 was: 

dp.l = 4.19590e + 00 dp. 5 = 6.84997e - 01 

dpJ2 = 3.03156e + 00 dp.6 = -1.69662e + 00 

dp. 3 = 4.25392e - 01 dp. 7 = 8.74300e - 02 

dp A = 8.41695e - 01 dp. 8 = -3.87025e + 00 

The largest percentage variation between these two optimal points is 65% for 
dp_7, and the largest variation is 0.0580 for dpJ2. These solutions are close one 
to another. 

The optimal point achieved when the optimization started in the trivial point 
and used epsilonF = 0 was: 

dp. 1 = 8.00974e - 01 dp. 5 = 4.72042e - 01 

dpJ2 = 1.34126e + 00 dp. 6 = — 3.33092e + 00 

dp. 3 = 1.30245e -I- 00 dp. 7 = 2.20404e - 01 

dp A = 5.26028e - 01 dp. 8 = -1.45150e + 00 
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The optimal point achieved when the optimization started in the trivial point 
and used epsilonF = 0.0002 was: 

dpA = 9.57795e - 01 dp. 5 = 7.17266e - 01 

dpJ2 = 9.82241e - 01 dpJo = -1.70196e + 00 

dp- 3 = 8.40876e - 01 dp. 7 = 4.88951e - 02 

dp A = 5.56358e - 01 dp. 8 = -1.07892e + 00 

The largest percentage variation between these two optimal points is 350% 
for dp.7, and the largest variation is 1.6290 for dp_6. 

It can be seen in Table 2.1 that for the F-14 problem the results were better 
when epsilonF=0 was used. The performance is presented in Figures 2.6-2.10, 

In Figure 2.7 it can be seen that the objective function, pitch bandwidth and 
time delay were improved at the expense of the tail rate soft constraint, which 
was pushed to its border. The speed of the airplane response cannot increase 
more, because this would violate the limit for the tail rate. Note that the YF-17 
spec is check only. Thus, CONDUIT ignores the fact that the phase portion of 
the spec is at Level 3. 

In Figure 2.8 it can be seen that the objective was improved at the expense 
of the attitude tracking soft constraint (Alpha error), which was pushed to its 
limit. The previous comment on the YF-17 spec applies here as well. 

In Figure ?? it can be seen that at the optimal point all objectives and con- 
straints are Level 1 excepting two of them (Alpha error and Gain Requirement) 
which exceed insignificantly the Level 1/2 border. Notice particularly that the 
theoretical model-following spec is satisfied. It is this spec that was used, not 
the YF-17, which is check only. 
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Figure 2.7: Grumman F-14, nominal point improved by CONDUIT, epsilonF=0 
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Figure 2.8: Grumman F-14, nominal point improved by CONDUIT, epsi 
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Figure 2.9: Grumman F-14, trivial starting point 
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Chapter 3 


Antenna Problem 

3.1 Description 

Usually a controller for an antenna is tested for tracking errors when following 
a command and for servo errors due to disturbances. The tracking and distur- 
bance rejection properties have to be traded off. The feedback loop is effective 
in suppressing the perturbation, but has weaker tracking properties. On the 
other hand the feedforward controller has good tracking properties, confirmed 
by both simulations and field measurements, but its ability to compensate the 
perturbations is insufficient. When comparing the properties of feedback and 
feedforward, one can conclude that by combining the two it is possible to im- 
prove both tracking and disturbance rejection properties. 

This problem is derived from a master’s thesis example, [5]. The problem is 
to design a servomechanism for controlling the line-of-sight (LOS) of an airborne 
electro-optical device about a particular axis of rotation. The controller’s most 
critical function is to isolate the LOS from perturbations induced by aircraft 
angular motions. So, only the feedback loop that assures perturbation rejection 
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will be discussed in the following. Aircraft angular motions produce disturbance 
torques which must be nullified by the servomechanism’s torque motor. The 
torque motor is driven from feedback elements mounted on the controlled device. 
These elements, a rate-measuring gyroscope, and an angular accelerometer are 
inertial in character and hence measure LOS in an absolute sense. 

Unlike the F-14 problem, the concern of this design is not the precision with 
which the system responds to command inputs. Rather, the goal here is to 
design a system whose response to disturbance inputs is minimized. The reader 
may find other instructive contrasts to the F-14 problem: 

1. As is typical in the servomechanisms field, where design via the time do- 
main methods of modern control theory has never really taken hold, the 
design criteria are specified completely in the frequency domain. 

2. Because higher frequencies are of interest here, more so than in the F-14 
problem, the design must consider dynamic characteristics of the actuator, 
mechanical structure, and feedback transducers themselves. These become 
part of the system ’’plant”. These issues are also important in aircraft 
SC AS and would be included in a more realistic SCAS design problem. 

3. Because of the above and the corresponding increase in the complexity 
of the specified form of the controller, the state dimension of the control 
system is significantly higher. 

4. Because the magnitude of the design parameters covers a much larger 
range, the problem is normalized by use of the ’’variation” specification, a 
feature of C-0 [1, pi. 6] that is not presently implemented in CONDUIT. 
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3.2 Block Diagram 


Figure 3.1 is the block diagram of an antenna control system. 

The components of the model are: 

• antenna dynamics 

• power amplifier and torque motor 

• angular accelerometer and accelerometer filter 

• rate gyroscope and gyroscope filter 

• control system, consisting of a lead-integrator, quadratic low-pass filter 
and quadratic lead-lag for the rate loop and lead-integrator and quadratic 
lag-lead for the acceleration loop 

The fixed parameter values of the plant shown in Figure 3.1 are: 
km = 0.801 2^ power amplificator and motor gain 
fm = 796 Hz motor electrical bandwidth 
J = 0.11 oz*in*s 2 gimbal inertia 
fr = 160 Hz gimbal resonance frequency 
far = 122 Hz gimbal anti-resonance frequency 
Zr = 0.015 gimbal structural zeta 
asf = 0.291 7J3- accelerometer scale factor 

IT 

fabw = 350 Hz accelerometer bandwidth 
faf = 500 Hz accelerometer filter bandwidth 
rsf = 2.263 7^ gyro scale factor 

a 

frbw = 75 Hz gyro natural frequency 
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LOS motion [miaofadians/(oz‘in)] 








zrt = 0.7 gyro zeta 

grfb = 10 ^ gyro filter gain 

frf = 482 Hz gyro filter bandwidth 

There are eleven parameters that have to be chosen to meet the desired spec- 
ifications: 

kl = rate-loop integral gain 

fl = negative of rate-loop lead pole 

f2 rate-loop second-order filter bandwidth 

f3 = rate-loop second-order lead-lag filter lag frequency 

f4 = negative of the acceleration-loop lead pole 

f5 = acceleration-loop second-order lag-lead lag frequency 

zl = rate-loop second-order lead-lag filter C 

z2 = acceleration-loop second-order lag £ 

al = rate-loop second-order lead-lap filter 

a2 = acceleration-loop second-order lag-lead filter 

In order to simplify the block diagram the design parameters used by CON- 
DUIT are formed from the previously listed parameters as indicated below: 
dpp. 1 = k 1 

dppJl = 2 * pi * f 1 
dpp. 3 = 2 * pi * fl 
dppA = 2 * pi * f 3 
dpp .3 = 2 * zl 
dpp. 6 = al 2 
dpp. 7 = k2 
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dpp - 8 = 2 * pi * f 4 
dpp-9 = 2 * pi * /5 
dpp_10 = 2 * z2 
dpp_ll = a2 2 

For this problem the design parameters are required to be positive. This insures 
stability and minimum phase for the lead-lag compensator in the rate loop and 
for the lag-lead compensator in the acceleration loop. 


3.3 Performance Objectives and Constraints 

The second part of the formulation of the control problem is the selection of 
the specifications. The specifications are the same as those in [5] except that 
two specs have been added. The specs have been rewritten in the appropriate 
form for CONDUIT. All the speciications are summarized below. They are all 
illustrated in Fig. 3.2. 

Objectives 

1. The objective is to minimize the motion of the LOS in response to a distur- 
bance torque input. This minimization is to be achieved over a disturbance 
frequency range from 0.1 Hz to 1000 Hz. The ’’good” and ’’bad” values for 
this quantity are 60 and 120 microradians/oz * in, respectively. 

torqspc code: 

closed loop, input 1 disturbance torque, output 1 LOS motion 
objective, u € [10 -1 , 10 3 ] Hz 

minimize, Level 1/2 = [-60 60][**gaf“], Level 2/3 = [-80 120][«*™tf“] 
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Constraints 


1. Stability criterion of the total feedback system (i.e., with both feedback 
paths active) is “all eigenvalues negative”. 

EigLcGl code: 

closed loop, from any input to any output 

[a,b,c,d]=closed-loop model, unstab equals sum of positive eig(a) if there 
are some positive eigenvalues or unstab equals max of negative eig(a) if all 
eigenvalues are negative 
hard constraint 

minimize unstab : Level 1/2 = 0.001, Level 2/3 = 0.002 

2. Stability margins of the feedback system are: 
gain margin > 6 [db] 

phase margin > 45 [deg] 

StbMgGl code: 

open loop, input 3 broken-loop input, output 4 broken-loop output 
hard constraint, u € [0.1,500] Hz 
maximize gain : Level 1/2 = 6 [dB], Level 2/3 = 4 [dB] 
maximize phase: Level 1/2 = 45 [deg], Level 2/3 = 35 [deg] 

3. Stability margin of the rate feedback loop must be insured by limiting the 
upper bound of the closed-loop gain. 

RaGMspc code: 

closed loop, input 2 gyro noise input, output 2 closed-loop rate gain 

soft constraint, u> e [10°, 10 3 ] Hz 

minimize, Level 1/2 = [0 1.4], Level 2/3 = [-0.5 1.8] 
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4. Stability margin of the acceleration feedback loop must be insured by lim- 
iting the upper bound of the closed-loop gain. 

AcGMspc code: 

closed loop, input 1 disturbance torque, output 3 closed-loop acceleration 
gain 

soft constraint, u € [10, 10 4 ] Hz 

minimize, Level 1/2 = [0 1.4], Level 2/3 = [-0.5 1.8] 

5. Two constraints are also included which are related to the simulation and 
not system performance itself. These involve the form of the two second- 
order compensators, one within the rate loop and the other within the 
acceleration loop. 

The intent of the first such constraint is to insure a lead-lag form for the 
rate loop compensator. This requires the constraint that al > 1. 
RComspc code: 
soft constraint 

maximize dpp . 6 = al 2 , Level 1/2 = 1.01, Level 2/3 = 1 

6. The intent of the acceleration loop constraint is to compensate for the 
structural anti-resonant-resonant characteristics of the gimbal. As such, 
the second-order filter must be of a lag-lead form which is insured by 
applying the constraint dL >1. For optimization it is better to use a 
smooth performance criterion. After several trials it was determined that, 
instead of using ^ > 1, it is better to use 0 < a2 2 < 1, or 0 < dppAl < 1. 
Denoting dp.ll as dpp.ll, the parameter is maintained positive. The 
upper limit is required by this spec. 
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Acmlspc code: 

soft constraint 

minimize dpp.ll, Level 1/2 = 0.99, Level 2/3=1 

3.4 Setup of the Initialization File 

The performance chart of this problem does not contain time analysis, only 
frequency analysis. Thus linear simulation will not be performed and there is no 
need to specify dt. This problem is too fast (the fastest mode is 796 Hz for the 
nominal starting point) to be practically simulated in time. 

epsilonF was set at 0.0002 (the default value) and 0 (machine precision), and 
tests were made for both values. 


3.5 Optimization 

For the antenna problem two starting points were considered. An acceptable 
starting point, called the nominal point, and an unacceptable starting point, 
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called the trivial point. The nominal point is defined by: 


kl 

~ 

>\ > 

00 

o 

t-H 

n 

= 

8.6 Hz 

f 2 

— 

417 Hz 

/3 

= 

250 Hz 

zl 

= 

0.707 Hz 

a 1 

= 

2 

k2 

— 

500 ^ 

/4 

= 

796 Hz 

/5 

= 

120 Hz 

z2 

= 

0.1 

a 2 

— 

0.817 Hz 


or equivalently approximately: 


dpp. 1 = 1.08000e + 02 

dpp. 7 = 5.00000e + 02 

dpp.2 = 5.40300e + 01 

dpp. 8 = 5.00132e + 03 

dpp. 3 = 2.62000e + 03 

dpp .9 = 7.54000e + 02 

dpp A — 1.57100e + 03 

dpp_10 = 2.00000e - 01 

dpp. 5 = 1.41370e + 00 

dpp .11 = 6.67500e — 01 

dpp.6 = 4.00000e + 00 
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and the trivial point by: 


kl = 50 - 

V 

/I = 5.0 Hz 
/2 = 100 Hz 
/3 = 100 Hz 
zl = 0.50 Hz 
a\ = 1.5 
k2 = 200 * 

V 

/4 = 200 Hz 
/5 = \00Hz 
z2 = 0.50 
a2 = 0.667 Hz 

or equivalently approximately: 


dppA = 5.00000e + 01 
dppJ2 = 3.14150e + 01 
dpp.3 = 6.28300e + 02 
dpp.4 = 6.28300e + 02 
dpp.5 = l.OOOOOe + 00 
dpp . 6 = 2.25000e + 00 


dpp.7 = 2.00000e + 02 
dpp.3 = 1.25660e + 03 
dpp£ = 6.28300e + 02 
dpp_10 = l.OOOOOe -I- 00 
dppAl = 4.44900e - 01 


The optimization results are in Table 3.1. The optimal point achieved when the 
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Table 3.1: Comparison between the optimal points obtained by CONDUIT for 
the antenna problem 



r 

Starting point 

Optimal using 
epsilonF=0 

Optimal using 
epsilonF =0.0002 

Nominal 

MAX. COST SOFT 

MAX. COST ^OFT 

MAX .COST JSOFT 


0.201512 

0.0132284 

0.0138598 

Trivial 

MAXJ1ARD 

MAX. COST JSOFT 

MAX.COSTJSOFT 


1.42554 

2.26952 

1.9644 


optimization started in the nominal point and used epsilonF — 0 was: 


dpp. 1 = 1.08002e + 02 
dpp.2 = 5.40251e + 01 
dpp. 3 = 2.62000e 4- 03 
dpp. 4 = 1.57100e + 03 
dpp .5 = 1.32943e 4- 00 
dpp Jo = 3.98639e + 00 

The optimal point achieved when the 

and used epsilonF = 0.0002 was: 

dpp. 1 = 1.08027e 4- 02 
dpp. 2 = 5.39746e 4- 01 
dpp _3 = 2.62000e 4- 03 
dpp A = 1.57100e + 03 
dpp _5 = 1.30668e -f 00 
dpp.6 = 4.06439e 4- 00 


dpp . 7 = 5.00002e + 02 
dpp . 8 = 5.00132e 4- 03 
dpp.9 = 7.54000e + 02 
dpp . 10 = 6.23329e - 02 
dpp.ll = 6.20870e - 01 

optimization started in the nominal point 


dpp. 7 = 5.00014e + 02 
dpp . 8 = 5.00132e + 03 
dpp.9 = 7.54001e + 02 
dpp . 10 = 4.88735e - 02 
dpp.ll = 6.83348e — 01 
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The largest percentage variation between these two optimal points is 21% for 
dp-10, and the largest variation is 0.0780 for dp.6. 

The optimal point achieved when the optimization started in the trivial point 
and used epsilonF = 0 was: 

dpp.l = 5.02402e + 01 dpp.7 = 2.00232e + 02 

dppJ2 = 3.07500e + 01 dppJi = 1.25660e + 03 

dpp- 3 = 6.28322e + 02 dpp. 9 = 6.28299e + 02 

dpp- 4 = 6.28306e + 02 dpp. 10 = 4.65478e - 01 

dpp- 5 = 6.93597e - 01 dpp. 11 = 1.01263e + 00 

dpp. 6 = 5.25785e + 00 

The optimal point achieved when the optimization started in the trivial point 
and used epsilonF = 0.0002 was: 

dpp- 1 = 5. 18130e + 01 dpp.7 = 2.02086e + 02 

dppj. = 2.57491e + 01 dppj& = 1.25658e + 03 

dpp- 3 = 6.28449e + 02 dpp_9 = 6.28298e + 02 

dpp A = 6.28578e + 02 dpp.10 = 7.9041 7e - 01 

dpp- 5 = 5.85324e - 01 dpp.l 1 = 1.00956e + 00 

dpp. 6 = 4.92225e + 00 

The largest percentage variation between these two optimal points is 69% for 
dp_10, and the largest variation is 5 for dp_2. 

It can be seen in Table 3.1 that for this problem the performance is better for 
the optimization using epsilonF=0.0002. 

The performance for the problem starting in the nominal point and optimized 
with epsilonF=0.0002 is presented in Figure 3.2 and Figure 3.3. By inspecting 
Figure 3.2 it can be seen that there is only one spec, AcGMspc, which does not 
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Figure 3.2: Antenna problem, nominal starting point 
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satisfy Level 1 at the starting point and it is improved through optimization, 
Figure 3.3. 

The performance for the problem starting in the trivial point and optimized 
with epsilonF=0.0002 is presented in the Figure 3.4 and Figure 3.5. By inspect- 
ing Figure 3.4 and Figure 3.5 it can be seen that the stability margins improved 
from the Level 3 (red zone) to the Level 1 (blue zone). The improved point has 
gain margin=6 and phase margin=45.5, being near the border of Level 1/2. The 
gain in stability was obtained at the expense of a soft constraint, Acmlspc, that 
couldn’t be satisfied. Optimization changed the unacceptable initial point to 
an acceptable optimal point. The unsatisfied spec, Acmlspc, is not important, 
because it is an artificial spec that requires the compensator in the acceleration 
loop to be lag-lead. Optimization showed that for these parameters a lead-lag 
works better. 


3.6 Optimization Attempts 

As a cautionary note, one has to be very careful with the constraints. 

Originally the antenna problem didn’t contain the EigLcGl and StbMgGl 
specs, and was optimized using a simplified performance chart containing only 
five specifications: RaGMspc, AcGMspc, RComspc.m, Acmlspc, torqspc. Fig- 
ure 3.6 and Figure 3.7 show the results of this optimization attempt. It can 
be seen that even though the system is stable at the starting point, it is driven 
to instability. This is obviously unsatisfactory, regardless of the values of the 
specifications actually included in the problem formulation. 

In addition to introducing the stability specs, another improvement done for 
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Figure 3.7: Antenna problem, nominal point is driven unstable when stability 
specs are check only, epsilonF=0.0002 
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this problem was the restriction that all design parameters be positive. Oth- 
erwise they can become negative, resulting in unstable or nonminimum phase 
compensators. 
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Chapter 4 


Sensitivity 


CONDUIT provides an optimal set of design parameters, 9, with respect to 
the given specifications. But the remaining question is if the controller model 
proposed by the designer is good enough. In a controller design problem, the 
sensitivity issue occurs in several ways: 

1. performance sensitivity to plant parameter changes 

2. performance sensitivity to controller parameter changes 

3. parameter sensitivity to small changes in the constraints 

An analogous sensitivity question arises in maximum likelihood estimation 
(MLE). Effective tools for analyzing sensitivity in the MLE situation are well 
known [16]. These tools have been implemented in a software package for aircraft 
identification and extensively tested and used [23]. 

The main focus of the work described here was on the extension of the MLE 
ideas and methods to CONDUIT and the control design problem. 
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4.1 The Different Forms of Sensitivity 

4.1.1 P erformance sensitivity to plant parameter changes 

Changes in the plant parameters are very much like disturbance signals and 
are generally addressed by insisting on robust controllers. Thus, the sensitivity 
of the design to changes in plant parameters is covered by the specifications, 
especially those relating to stability margins. 

4.1.2 Performance sensitivity to controller 
parameter changes 

Min max optimization is trying to minimize the maximum objective (constraint) 
at every iteration. This function is called MAX-ACTIVE in C-O. The MAX-ACTIVE 
function, called MAX-HARD (see section 1.2) in phase 1, MAX_COST_SOFT 
in phase 2, or MAX.COST in phase 3, is the function that is actually being 
minimized. The possible situations when optimization stops are presented in 
Figures 4. 1-4.8. Note that, in these figures, specifications are indicated by thin 
lines while the active spec (MAX-ACTIVE) is indicated by a heavy line. 

• One or some hard constraints cannot be satisfied, and the optimization 
stops in Phase 1, Figure 4.1 and Figure 4.2. For these cases, the sensitivity 
computation is unnecessary, because a feasible set of parameters and an 
optimal point do not exist. 

• All hard constraints are satisfied, one or more soft constraints cannot be 
satisfied, and the optimization stops in Phase 2, Figures 4.3-4. 5. 
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Figure 4.2: The minimum is Phase 1 
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It is also possible that the optimization stops in Phase 2 when an objective 
competes with a soft constraint. But in this situation the optimization 
problem is not well formulated and it is better to change that objective so 
it is satisfied. This way the objective will allow the soft constraint to be 

improved first, and then, in Phase 3, the objective will be improved to its 
best value. 


In Figure 4.5 it can be seen that in the region where the 


hard constraint >1, (4 1 ) 

MAX_ACTIVE equals the hard constraint even if there is a larger soft 
constraint. 


When the minimum is in Phase 2 the optimization is considered acceptable. 

• When all hard and soft constraints are met, the optimization stops in Phase 
3, as shown in Figures 4.6, 4.7, and 4.8. 

Figure 4.8 shows a situation where the constraints are met (< 1), in the 
region surrounding 0; MAX-ACTIVE equals the objective even if there is a 
constraint with a larger value. This case, in which the optimization pushes 
a constraint to its border by minimizing an objective is the most common 
case. 

When the minimum is in Phase 3 the optimization is considered satisfac- 
tory. 
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Figure 4.6: The minimum is Phase 3 
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Figure 4.7: The minimum is Phase 3 
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Figure 4.8: The minimum is Phase 3 
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These minimum points can be characterized by considering the local geome- 
try in two classes: 

A) When the minimum occurs at a point similar to those shown in Figures 
4.3 and 4.6, it is equivalent to the solution to an unconstrained optimization 
problem. MAX-ACTIVE has: 

• Gradient zero on all design parameter axes at the optimal point 

• Hessian positive semidefinite at the optimal point 

and J, the function used in sensitivity computation, equals MAX-ACTIVE. 

B) When the minimum looks as in Figures 4.4, 4.5, 4.7, and 4.8 
MAX-ACTIVE has: 

• Gradient discontinuous on at least one design parameter axis 

• Hessian undefined. 

The sensitivity test implementation includes the gradient computation of all 
optimization specs and of J, at the optimal point. In case the gradient of J 
is zero, then the optimal point is case A) and insensitivities, correlations, and 
Cramer-Rao bounds can be computed as in section 4.2. In case the gradient of 
J is discontinuous, then section 4.3 applies. 

4.1.3 Parameter sensitivity to small changes in the con- 
straints 

The sensitivity of the optimal solution to small modifications in the constraints 
is an important aspect of practical applications. The sensitivity theorem will be 
presented in section 4.4. 
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4.2 Sensitivity Tools for Use When the Gradi- 
ent of MAX_ACTIVE is zero at the optimal 
point 

4.2.1 Insensitivity Ellipsoid 

• Estimation 

In estimation theory, the confidence region of an estimator is often evalu- 
ated using the probability density function of the estimate, p§{0), i.e. the 
probability that the estimated 0 is actually 0. If the maximum likelihood 
estimation technique is used to estimate a parameter vector 0, this prob- 
ability density function is approximately Gaussian, unbiased, and with 
covariance matrix equal to M(0t)~ 1 '- 

p § (0) 2 [(27r) n |M(0 t ) -1 1 l' 1 / 2 exp {-±(0 - 0 t ) T M(6 t )(6 - 0 t )} (4.2) 

where M is the Fischer information matrix: 

M(0 t ) = £{[V fi( \np(Z\0 t )}{Vl lnp(Z|0 f )|0 t ]} = ~E{V 2 0t \np(Z\0 t )\0 t } 

(4.3) 

where Z represents the observations and 0 t denotes the true value of 0 [17]. 
The equation 

(0 - 0 t ) T M(0 t )(0 - 0t) = Constant (4.4) 

describes isometric surfaces (ellipsoids) centered at 0 t , with the same value 
for the probability density. For Constant equal one, the volume integral 
over the interior of the ellipsoid is 63.2% and is equal to the probability 
that 0 lies within the ellipsoid. From the practical point of view 0 t is not 
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available and the above ellipsoid cannot be drawn. This problem is handled 
by reversing the roles of 9 t and 0, because M(9) can approximate M(6 t ). 
The ellipsoid defined by 

(9-9) T M(0)(0-9) = 1 (4.5) 

is called the uncertainty ellipsoid. For maximum likelihood estimation, 
if the negative log likelihood function J is: 

J ( d ) = \ Efrft) ~ z{U)] T R- l [zo{ti) - z(ti)] (4.6) 

the information matrix can be closely approximated by 

N 

H = 'H'V 9 zl(t i )R- l Vlz 9 (t i ) (4.7) 

1=1 

z is the observation vector, z is the predicted value of 2 if 0=9 , R is the 
covariance matrix, N is the number of measurements of 2 . The H-matrix 
is identical to the dominant term of the Hessian matrix VjjJ(0). 

• Optimization 

The estimation problem is to find the model parameters, 0 , that minimize 
the negative log likelihood function J. This function expresses the difference 
between the predicted output based on 0 and the measured output. 

The optimization problem is to find the controller parameters, 0, that 
minimize the MAX_ACTIVE function, denoted also by J. This function 
expresses the difference between the predicted performance based on 9 
and the desired performance. MAX_ACTIVE is called MAX. HARD in 
Phase 1, MAX.COSTJSOFT in Phase 2, or MAX.COST in Phase 3 of the 
optimization. Because there is no probability in the optimization case 

M(0) = VjJ(9). (4.8) 
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For case A) MAX-ACTIVE can be expanded about 9 t as follows: 

J{6) = J{0 t ) + \(0~ 0t) T V 2 e J(6 t )(6 - 9 t ) + 0{(6 - d t f) (4.9) 

9 t represents the parameters that achieve the best performance and 9 the 
optimal parameters as determined computationally. Suppose that J{9)- 
J(9 t ) is less than some constant. This situation is referred to as a sub- 
optimal solution. This can be sufficient where the exact optimum is not 
required, but performance is required to be within a specified range of the 
optimum. If higher order terms are ignored, this condition becomes 
1(0 - 0 t ) T VlJ(9 t )(9 - 9 t ) < Constant and can be approximately as before 
with 

~(0 t _ 0) T V 2 6 J{9)(9 t -0)< Constant (4.10) 

2 

This ellipsoid is essentially the same as the uncertainty ellipsoid if the con- 
stant is chosen to be But because it is described using MAX-ACTIVE 
and not the negative log likelihood, we named it the insensitivity ellip- 
soid. 

In the CONDUIT optimization problem the fit error, S7 2 e J{9), often cannot 
be approximated by H. If MAX-ACTIVE were computed as the sum of er- 
rors of performance functions then the Hessian of MAX-ACTIVE could be 
approximated by the sum of square gradients of all the optimization specs. 
But MAX-ACTIVE is computed as the max error of performance func- 
tions and shows the behavior of only one performance function at a time, 
without any influence of the others. Thus the Hessian of MAX-ACTIVE 
generally shows the curvature of only one performance function at a time. 
Using the squared gradient of MAX-ACTIVE to approximate the Hessian 
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frequently results in a singular matrix so the estimation theory approxi- 
mation cannot be applied. These limitations do not apply if an objective 
is formed by summing a collection of specs. It is recommended that the 
summed objective be used in CONDUIT. 

Moreover, if M is approximated with the sum of squared gradients or the 
sum of Hessians of all optimization specs, then there is not certitude that 
M is positive semidefinite. The specs are randomly convex, concave or 
neither at the optimal point. But the Hessian of J is positive semidefinite, 
as will be shown in section 4.2.5. 

The insensitivity ellipsoid can be thought of as the result of minimizing 
the MAX_ACTIVE. This minimization is a form of a sensitivity analysis 
based on the principle that two values of 0 cannot be reliably distinguished 
unless they result in sufficient differences in the MAX-ACTIVE function. 

Figure 4.9 shows the insensitivity ellipsoid in the case of two unknowns. 
The interior of the insensitivity ellipse is shaded. Most of the concepts 
discussed below can be illustrated by using such a two-dimensional ellipse. 

Optimization problems solved with CONDUIT have typically more than 
10 parameters; that is, the ellipsoid is more than 10 dimensional. Because 
it is hard to visualize the information on the ellipsoid, there are three other 
formats for displaying it: insensitivity, correlation, and Cramer-Rao bound 

[23], 
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Figure 4.9: Two-dimensional insensitivity ellipsoid 
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4.2.2 Insensitivity 

Insensitivity is a measure of how much one parameter can be changed from the 
achieved optimal value without causing the performance to degrade by more than 
a given amount. To simplify the computation, J is approximated by its Hessian 
(second order approximation), M, at the optimal point. This way, instead of 
recomputing J for parameter changes, only two more evaluations in addition to 
the optimal point evaluation are needed. 

An expression for the insensitivity can be obtained from the equation for the 
insensitivity ellipsoid, Equation 4.10. Investigate changes in 8 along the direction 
of the ith parameter, that is 8 — 6 = ae if where ej is a unit vector in the direction 
of the ith parameter and a is an arbitrary scalar. To stay inside the ellipsoid, 
the following must be true: 

a 2 e[Mei< 1 (4.H) 

Since e t Me x > 0 because M is positive semidefinite (see section 4.2.5), 

then 

M < ( Mu )~ l/ 2 (4.12) 

The maximum magnitude of a is given by the equality, which occurs at the inter- 
section of the ej axis with the ellipsoid. This value is defined as the insensitivity 
with respect to 0 t . Figure 4.10 illustrates the geometric interpretation of the 
insensitivity for a problem with two unknowns. For convenience, the origin has 
been redefined to be the center of the ellipse, which is the maximum likelihood 
estimate. The insensitivities of 6 X and d 2 are labeled as / x and / 2 . 

Remark 1. The insensitivity is a reasonable measure of accuracy only when a 
single parameter is being estimated, because it ignores any effect of correlation 
between parameters. When several parameters are estimated, the error band is 
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Figure 4.10: Geometric interpretation of insensitivities and Cramer- Rao bounds 
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a result of both insensitivity and correlations. Insensitivity gives only the lower 
bound for the error band, but correlation effects tend to increase the error band 
so much that the insensitivity is virtually useless as an indicator of accuracy. 

Remark 2. In order to compare insensitivities of different parameters to each 
other, they must be scaled with the parameter values !*■. 

4.2.3 Cramer- Rao Bound 

The Cramer- Rao bound is also based on the insensitivity ellipsoid [17], The 
standard deviation of the ith parameter is For maximum likelihood 

estimates, the Cramer-Rao inequality gives the approximation for the co- 
variance matrix and the quantity V /(7F T )~ is called the Cramer-Rao bound. In 
this thesis J is different from maximum likelihood estimation and the Cramer- 
Rao bound will be defined to be equal to the standard deviation. The Cramer- 
Rao bounds and insensitivities are closely related. While the insensitivity is 
the conditional standard deviation of the parameter estimate, given that all the 
other parameters are known, the Cramer-Rao bound is the unconditional stan- 
dard deviation. The insensitivity is ^/(M*)- 1 and the Cramer-Rao bound is 
The geometric relationship between these quantities is the following. 
The insensitivity is the solution to the following optimization problem: maximize 
\9i-Q x \ while staying within the insensitivity ellipsoid and holding all of the other 
parameters fixed. The Cramer-Rao bound is the solution to the same problem 
without holding the other parameters fixed. The removal of the restrictions on 
the optimization is directly analogous to removing the statistical conditioning. 

A two-dimensional plot of the geometric relationship between the two measures 
is shown in Figure 4.10. The insensitivities of 0 X and 6 2 are labeled I Y and / 2 , 
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while the Cramer- Rao bounds are labeled CR\ and Ci? 2 - 

The proof of the affirmation that the Cramer-Rao bound is the solution to 
the optimization problem above comes from the next theorem. 

Theorem 4.1 Given a fixed vector y and a positive definite symmetric matrix 
M, then the maximum of x T y, subject to the constraint that x T Mx < 1, is given 



Apply this theorem to the insensitivity ellipsoid, where y is e,, a unit vector 
along the axis of the ith parameter, and x is 0 — 0. The maximum of 0, — 0, = 
(0 — 9) T e{ is 

M -1 ej = (4.13) 

and this is CR, , the Cramer-Rao bound. CRi is the largest projection on the 0, 
axis of the vectors with origin in 0 and arrow on the ellipsoid. It is caused by 
the vector : 

9-9 = M ~' C| (4.14) 

This projection shows the largest error that can occur for 0, without affecting per- 
formance. This largest error for 0, is unperceived only if all the other parameters 

A ■ r -m~- J 

have particular values, so that together they form the vector 0 — 0 = —j==&=. 

Remark 1. In order to compare Cramer-Rao bounds of different parameters 
to each other, they must be scaled with the parameter values, 

Remark 2. To see the correlations that may cause a big Cramer-Rao bound, 
CRi , the vector 0-0 that has this projection on 0* must be decomposed in the 
vector space formed by all insensitivities vectors: J x * e x , e n , where e, is 

a unit vector along the axis of the ith parameter. This is equivalent to have an 
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ellipsoid that cuts all axes at 1, all insensitivities being 1. Figure 4.11 shows the 
projections p x , .. , p n . 


8-8 = 


M 


-i 


ei 




= \Pi P 2 ■■■ p n ] T = [ai * h a 2 * /; 


2 ■■ .a, 


* Inf ( 4 . 15 ) 


so that 


[ a l °2 • • • a-nf = = 


CRi 


( 4 . 16 ) 


where M J (:, i ) denotes the ith column of M -1 and T = / 2 , . . . /„]. The 

large components of Q CR% show the parameters that are correlated with Si. In 
this case the table showing the correlations is not symmetrical. 


We can compare the projections of the vector 8-8 in the vector space 
CR X *e u ... CRn * e n , where e { is a unit vector along the axis of the ith pa- 
rameter. This is equivalent to having an ellipsoid that can be enclosed in an 
n-dimensional hypercube with all sides 2, all CR , being 1. Then the correlation 
table is symmetrical, with all diagonal elements equal to 1. Both scalings (with 
insensitivities and with Cramer-Rao bounds) are equally good, and the first one 
is presented in the example in subsection 4.2.8. 


4.2.4 Eigenvalues and Eigenvectors 

The parameter correlation can be seen also by studying the eigenvalues and 
eigenvectors of M, normalized to have unity diagonal elements. The unnormal- 
ized M matrix is of very little use in studying correlation, because scaling effects 
tend to dominate. The normalization is done by dividing every row i and column 
i by ^/M(i, i). The resulting ellipsoid will intersect every axis at 1, having the 
insensitivity of all parameters equal 1. 
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Figure 4.11: Projections 
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If Aj, A 2 , • • • A n are the eigenvalues of the normalized M matrix, then the 
diagonals of the ellipsoid are 2^-. Small eigenvalues of the normalized M matrix 
indicate correlations among the parameters with significant components in the 
corresponding eigenvectors. 

Eigenvalues and eigenvectors values of the Hessian of MAX_ACTIVE are 
interpreted in the context of parameter correlation in an example at the end of 
subsection 4.2.8. 

4.2.5 Positive Semidefiniteness 

In order to compute the Cramer-Rao bounds, a positive definite M is necessary. 
The question is if the Hessian of MAX-ACTIVE is always positive. The answer 
is usually yes, but under certain circumstances it could also be zero. 

As can be seen in Figure A), the Hessian is positive definite. It is also pos- 
sible that MAX-ACTIVE reached a minimum in some direction, but does not 
change at all in some other directions. 

If MAX-ACTIVE is completely insensitive in a direction that coincides with a 
parameter axis, then the insensitivity ellipsoid is degenerate on that parameter 
axis, Figure 4.12. Therefore the Hessian is singular, insensitivity is infinite for 
that parameter and the Cramer-Rao bounds are also infinite. This may happen 
because that parameter does not change locally, an acceptable behavior, or glob- 
ally, in which case it should be eliminated. 

If MAX-ACTIVE is completely insensitive in a direction that does not coin- 
cide with a parameter axis, then the insensitivity ellipsoid is degenerate in that 
direction, but still cuts the axes, Figure 4.13. Therefore the Hessian is also sin- 
gular, insensitivities are finite, and the Cramer-Rao bounds are infinite. This 
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Figure 4.12: Ellipsoid is degenerate 





shows correlation, and one of the parameters involved in the correlation should 
be eliminated. 


4.2.6 Hessian Formulae 


To compute the Hessian one can use, among others, the following three methods: 
• Method one (faster): 

Hess jk (f, h) « \ * | / (xj + h, x± ± h) - 2f{x^, x k ) + f{xj - - h) 

fixj d~ h., Xfc) 2/(xy,x^) 4- /(x^ h, x k ) 

/r 

f (,•%] i d” ^ ) 2/(Xj , Xfc) T f (Xj , Xfc h) 1 

1? J 

(4.17) 

For this method one needs in general N 2 4- N function evaluations, where 
N is the number of the parameters. First the diagonal terms are evaluated 
using 2 * N function evaluations and then the upper diagonal elements 
needs only N 2 — N function evaluations. 

The truncation error for off-diagonal elements is 

+ 3 a^/ a foj.Sfc) + 2 a,f /ji fxjpXfc) + 

[ dx*dxk oxjdxi, dxjdx k J ( 4>18 ) 

Higher order terms 

The truncation error for diagonal elements is 


h 2 

T2 


12 


a 4 / 


dx* 


4- Higher order terms 


The roundoff error for off-diagonal elements is 

12c// 

2h 2 

The roundoff error for diagonal elements is 

4c// 

h 2 


(4.19) 


(4.20) 


(4.21) 
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• Method two (slower): 


Hessjk(f, h) 


f “h x k d" h) f {xj + /i, — /i) 

j 

/ (x ? — h, Xfc + /i) — / (z, — h, xjt — A) 


4/r 


(4.22) 


For this method one needs in general 2 * N 2 function evaluations, where 
N is the number of the parameters. 


h 2 [ d*f d*f 

J [dxpr k ^ Xj,Xk) + d^d4 {Xj ' Xk) 


The truncation error for off-diagonal elements is 

-I- Higher order terms (4.23) 
The truncation error for diagonal elements is 

-I- Higher order terms (4.24) 


h? 

12 


a 4 / 


dxlj 


The roundoff error for off-diagonal elements is 


4e// 

4/i 2 

The roundoff error for diagonal elements is 

4 €ff_ 

h 2 


(4.25) 


(4.26) 


• Method three: 


Hess jk (f,hj,h k ) « 


/ ( x i hj , x k + /it) — f(xj + hjjXk ~ h k ) 

4hjh k 

£i ~ h u Xk /tfc) — f(xj — hj, x k — h k ) 
4hjh k 


F 


4.27) 


where hj = h* max(l, abs (xj)). 


For this method one needs in general 2 * N 2 function evaluations, where 
N is the number of the parameters. 
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The truncation error is 


h lJL Lfa. Xt) + *L«V_( j 

6 dxtdxu'' >' k) 6 dx,dxV ’ ' 


+ Higher order terms (4.28) 


h for the Hessian equations is computed in subsection 4.2.7. 

Comparing the first two methods, it can be seen that for h of the order 
of machine precision, they are almost equivalent. In the CONDUIT case, h is 
bigger because the performance accuracy is lower. A good choice for h will be 
seen in the next subsection. For a bigger h, the second method for computing 
the Hessian is more accurate than the first one. Both the truncation error 
and the roundoff error are smaller for the second case. The truncation error 
of method one includes the term i x ji x k)- This term is usually bigger 

than $ - ££ - k {xj,Xk) or d £*£ x 3 {xj,x k ) because it contains lower-order partials in 
any parameter. In case that the third order partials are zero, this term may 
still exist. The strength of the second method consists in its symmetry. Both 
the first and the second method were coded and tested. The second method is 
recommended and used in the example presented in subsection 4.2.8. 


4.2.7 Results on Accuracy of Hessian Computation 

There are two sources of error in Hessian equations: truncation error and round- 
off error [18]. 


The roundoff error: With h exactly representable, one has a roundoff 
error in method 2: 

4c/ ^ (4.29) 


cr ~ 4 h 2 

where tf is the fractional accuracy with which / is computed 
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• The truncation error: 


,2 r^4 




dyfa.g*) dVfc, 
dxjdxl dxjdxk J 


+ 0(/i 4 ) 


(4.30) 


To minimize e r + e t , e r — e t . Therefore the optimal choice for h is 



where x c « (/// ) 7 is the “curvature scale” of the function f, or “character- 
istic scale over which it changes. In the absence of any other information, one 
often assumes x c = x (except near x = 0). 

In the CONDUIT case h = (0.0002) (1/4 > = 0.1189. 

Remark 1. In order to have tj — 0.0002 precision for specs they have to be read 
from the Pcomb chart with at least 4 significant digits instead of 2, and we made 
this modification in pcomb.c. 

4.2.8 Example 

An example of the case where the gradient of MAX_ACTIVE is zero is the opti- 
mal point found for the F-14 problem, starting at the trivial point and optimizing 
with epsilonF = 0.0002. The Pcomb chart shows that the optimal point is the 
local minimum of the spec aerrspc, which does not compete with any other spec, 
as can be seen in Figure 4.14: 

The sensitivity analysis was done for this point, using the second method for 
the Hessian computation and step h=0.1189. The results found are displayed in 
Figures 4.15, 4.16, and 4.17. 
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Figure 4.14: Grumman F-14, performance display in the Pcomb chart 
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Figure 4.16: Grumman F-14, insensitivities 
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Figure 4.17: Grumman F-14, correlations 
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The largest Cramer-Rao bound corresponds to dp_7 and this is a consequence 
of a large insensitivity of the performance to this parameter. For the F-14 
problem this does not necessarily mean that dp_7 is unimportant, but there is 
no spec to drive the parameter in this phase of the CONDUIT design run. dp_7 
is the time constant in the stick prefilter and it is tuned to reject the noise added 
to the pilot command. Because in this simulation there is no noise added to 
the cockpit, performance is insensitive to dp_7. Sometimes a parameter appears 
insensitive because the solution is very robust in that parameter at the optimal 
point. 

A much more reliable result is the correlation result. Figure 4.17 shows the 
correlation of design parameters having Cramer-Rao bounds larger than 20%. 
Every shown column corresponds to such a parameter and shows the correlations 
of this parameter with all the other design parameters. For this particular case 
dp_3, dp_6, dp_7, and dp_8 have Cramer-Rao bounds larger than 20% and they 
are analyzed for correlation. As was explained, dp_7 has Cramer-Rao bound large 
as a consequence of its insensitivity. But the others are sensitive parameters and 
their large Cramer-Rao bounds are a consequence of their correlations. Figure 
4.17 shows that : 

- dp_3 is correlated to dp_8 

- dp_6 is correlated to dp_7 and dp_8 

- dp_8 is correlated to dp_6, dp_3 and dp_7. 

This suggests that the elimination of dp_8 may be a good approach. dp_8 is the 
integrator gain and dp_4 and dp_5 can do its job. 

For the same optimal point we computed the eigenvalues and the eigenvec- 
tors of the scaled Hessian of the MAX-ACTIVE. The eigenvalues are shown 
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Figure 4.18: Grumman F-14, Hessian’s eigenvalues in the optimal point 

in Figure 4.18 The corresponding eigenvectors are displayed in Table 4.1. The 
eigenvalues plot shows that the largest eigenvalues are the seventh and eight, 
corresponding to dp_7 and dp_8. In the eigenvectors Table 4.1, every column is 
an eigenvector corresponding to a dp and shows the correlation of the dp with 
all the others. Note that dp_7 is highly correlated to dp_6 and that dp_8 is cor- 
related to all dp’s but dp_7. These results are in accordance with the results of 
the Cramer-Rao ellipsoid. 
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Table 4.1: Grumman F-14, Hessian’s eigenvectors in the optimal point 


1 

vector 

2 

vector 

3 

vector 

4 

vector 

5 

vector 

6 

vector 

7 

vector 

8 

vector 

0.5401 

-0.2344 

0.1045 

0.6348 



0.1197 

-0.4522 

-0.3502 

0.1593 

-0.1088 

0.7643 



-0.1182 

0.4670 

-0.4322 

-0.5166 

0.1063 

0.0336 



-0.0098 

-0.3733 

0.3676 

0.0898 

0.7463 

-0.0468 



-0.1344 

0.3408 

0.4760 

-0.3160 

-0.5830 

-0.0851 

0.3545 


-0.3376 

0.2922 

-0.0445 

-0.2497 

0.0209 

0.0013 

0.2438 

-0.6692 

0.5790 

0.3045 

0.1874 

0.3423 

-0.2363 

-0.0159 

0.2613 

0.5206 

0.6723 

-0.0192 

-0.0020 

-0.6044 

0.1132 

-0.0454 

-0.4997 

0.4174 

0.2289 

0.3788 


4.3 Gradient of MAX_ACTIVE is discontinu- 
ous at the optimal point 

Consider the case of the constrained minimization problem 

minimize f(x) 

subject to h(x) = 0 (4.32) 

p(x) < 0 

f : R n —> R,h : R n : RJ 1 — ► R p , all C l , and the optimal point 

x*. In case that there are more objectives, f(x) is the maximum of them. If 
the optimization stopped in Phase 1 or Phase 2, when not all hard and soft 
constraints are satisfied, then the following discussion does not apply. It is true 
only for optimal points that satisfy all the constraints. 

If the optimal point x* is found in case B), the gradient of MAX-ACTIVE 
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is discontinuous, but the gradients of the objective and constraints satisfy the 
Kuhn-Tucker Conditions [21] [22]. 

Definition 4.1 [21] Let x* be a point satisfying the constraints 

h ( x *) = 0, g{x*) < 0 (4.33) 

and let J be the set of indices j for which gj{x m ) = 0. Then x* is said to be 
a regular point of the constraints 4.33 if the gradient vectors V/i t (x’), Vg^x'), 
1 < i < m, j £ J are linearly independent. 

Theorem 4.2 (Kuhn-Tucker Conditions) [21] Let x* be a local minimizer 
for the problem 4-32 and suppose x * is a regular point for the constraints. Then 
there is a vector X £ Rm and a vector p £ Rp with p>0 such that 


V/(x*) + XVh(x*) + pVg(x*) = 0 

(4.34) 

pg{x*) = 0 

(4.35) 

h(x*) = 0 

(4.36) 


The Kuhn-Tucker Conditions still hold when the regularity conditions for x* 
are replaced by a weaker condition, i.e. the Kuhn-Tucker constraint qualification 
(KTCQ) satisfied at x* [22]: 

coTC(x*,ft) = {h\ < Vgj{x*),h>< 0 Vj 6 J(x*)}(= S(x*)) 

where the tangent cone is defined by: 

TC(x,Q) = {h £ i? n | 3o(),Z > 0 s.t. x + th + o(t) £ Q 
Vf £ (0, t], ^ 0 as t — > 0, t > 0} 
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and J7 is defined by Equation 4.35 and Equation 4.36, co(X) is the union of all 
convex combination of finite subsets of X, and S(x*) is called the set of “first- 
order feasible “ direction [22]. 

By inspecting Equation 4.34 for the very simple and frequent case when at 
the optimal point the objective pushes one constraint to its border, 

V/(x*) + fijVgj(x*) = 0, j e J{x *) 

it can be seen that the scalar product of the two normalized gradients is pj = — 1 . 
In this case the other components of p, corresponding to non-active constraints, 
are zero. Analyzing the normalized gradients of the objectives and constraints 
at the optimal point, the active constraints can be recognized. 

In case that there are more active constraints at the optimal point, their gradients 
satisfy Equation 4.34, but it is harder to identify them among all the active 
and non-active gradients. When the analytical formulae for the objectives and 
constraints are known the solution of the system can be found by defining various 
combinations of active constrains, solving the system of n+p+m equations with 
n+p+m unknowns, and checking the signs of the resulting Lagrange multipliers. 
fj.j must be positive or zero. 

Second-Order Conditions exist also for problem 4.32 at the optimal point. 

Theorem 4.3 (Second-Order Necessary Conditions) [21] [22] Suppose 
the functions f, g, h G C 2 and that x* is a regular point of the constraints 4-33. 
If x* is a local minimizer for problem 4-32, then there is a X € R 171 , p € R p , 
p > 0 such that Equations 4-34 and 4-35 hold and such that 

d 2 L 

<t,—{x*,\,p)t>>0 (4.37) 


no 



on the tangent subspace of the active constraints at x* , with 

L(x, A ,p) = f (x) + A h(x*) + pg(x*) 

These conditions reveal the existence of a positive semidefinite Hessian on the 
tangent subspace: 

. ,dh , 

Vi €{ fc {x)t = 0 ' <Vg j (x*),t>=OVjeJ(x*)} ( 4 . 38 ) 

Testing if the Hessian of the Lagrangian L can be used instead of the Hessian 
of MAX-ACTIVE in computation of the Cramer-Rao bounds, insensitivities, 
and correlations, in a constrained optimal point, remains a suggestion for future 
work. 


4.4 Parameter sensitivity to small changes in 
constraints 

A sensitivity theorem exists to make sure that in the usual cases (functions are 
C 2 ) if an optimal solution satisfying all the constraints exists, then after small 
variations of the constraints an optimal solution continues to exist. 

Theorem 4.4 [21] [22] Let f,g,h E C 2 and consider the family of problems 

minimize f(x) 

subject to h(x) = c ( 4 . 39 ) 

g{x) < d 

Suppose that for c=0, d=0, there is a local minimizer x* that is a regular point 
and that, together with the associated Lagrange multipliers, \, p > 0, satisfies the 
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second-order sufficiency conditions for a strict local minimum. Assume further 
that no active inequality constraint is degenerate. Then for every (c, d) € R m+P 
in a region containing (0,0) there is a solution x(c,d), depending continuously 
on (c,d), such that x(0, 0) = x* , and such that x(c,d) is a local minimum point 
of 4-39. Furthermore, 

V c /(x( 0,0)) = -A (4.40) 

V d (x(0,0)) = — /i (4.41) 
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Part IV 

Conclusions and Suggestions for 

Future Work 
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In this report, two multivariable control system design problems were ana- 
lyzed and optimized using CONDUIT. This section will summarize the conclu- 
sions that are drawn from the results presented in Chapter III and then suggest 
possible paths of continued research. 

First we must emphasize that CONDUIT is an invaluable tool for aircraft 
control system design design, that allows a quick analysis and design of different 
control problems. Every change in the design can be executed instantaneously 
because there is no need for code writing. 

Although the two discussed problems are very different, they were satisfac- 
torily optimized. The system in the F-14 problem is a rather slow system and 
was analyzed especially in the time domain, while the system in the antenna 
problem is fast and was analyzed only in the frequency domain. 

Before starting the optimization process the models and the specifications had 
to be changed as follows: 

• F-14 problem: the block diagram was improved by adding a wind gust 
model and the performance chart was supplemented with three specs - 
Bandwidth & Time delay, RMS g’s at the pilot caused by turbulence, and 
the eigenvalues test 

• Antenna problem: stability was insured by adding an eigenvalues test spec 
and stability margins spec, and by requiring positive design parameters. 

The value of the discretization sample rate for the time simulation, the value 
of the gradient step-size, and different initialization points were tested. The 
sensitivity issues were addressed. The performance sensitivity to controller pa- 
rameters changes was detailed and parameter insensitivity and correlation were 
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derived as an approximation of the estimation theory. These were illustrated for 
an optimal point achieved for the F-14 problem. 

Some problems still remain open and should be addressed in future works: 

• How can changes in the sample rate affect the optimal design and how 
sensitive is the optimal design to them ? 

• Why does the optimal point vary slightly when starting the optimization 
in different points close to each other ? 

• Is it possible to use the gradients or the Lagrangian instead of the Hessian 
of J in the constrained optimal point to compute Cramer-Rao bounds, 
insensitivities, and correlations of the parameters ? 
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